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, Abstract 

o 



We consider the chemostat model with the substrate concentration as the single measurement. We 
^ ■ propose a control strategy that drives the system at a steady state maximizing the gas production without 

the knowledge of the specific growth rate. Our approach separates the extremum seeking problem from 
the feedback control problem such that each of the two subproblems can be solved with relatively simple 
algorithms. We are then free to choose any numerical optimization algorithm. We give a demonstration 
for two choices: one is based on slow-fast dynamics and numerical continuation, the other is a combination 
of golden-section and Newton iteration. The method copes with non-monotonic growth functions. 
^ Key-words, adaptive control, self-optimizing control, parameter optimization, biotechnology. 

^ ■ 1 Introduction 

■ The control design of chemostat models have been extensively studied in the literature, with the objective 
I to provide efficient and reliable control strategies for industrial applications, such as biotechnological or 

■ pharmaceutical processes [9l [TOl [15] . A common task is to drive a continuous stirred tank bioreactor to a 
set-point that optimizes a objective function, for instance the gas production rate [111 I43j . In this work, we 
consider the chemostat model with a single strain [33] 

- T— I 

S = -^J.{s)b + u{sin- s) 

^ : b = li{s)b-ub ^ > 

where the state variables s and b denote the substrate and biomass concentrations, respectively (in these 
equations the yield coefficient has be chosen equal to one without any loss of generality) . The concentration 
of substrate in the feed is denoted by Sin, and the dilution rate u > is the manipulated variable. The 
production rate to be maximized is given by 

r = /i(s)6 . 

In such bioprocesses, if often happens that the growth kinetics /i(-) is unknown (or poorly known) and possibly 
evolves slowly with time or changes of environment (temperature, pH...). Consequently, the robustness of 
the control strategies with respect to uncertainties on /i(-) is a crucial issue for real applications. Many works 
have been done for the on-line estimation of the growth function [3T1 [T31 [33 [3H1 [351 [Z] and the robust 
stabilization of such processes about a given reference point [T31 [HI [3S1 [SZl [30] , but there are comparatively 
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much less works that concern the maximization of an objective under model uncertainties [TB]. Such issues 
are typically addressed by the design techniques of extremum seeking controls [50l [H [28l [26l |3l IH |45l |46] . 
Roughly speaking, these techniques consist in adding an excitation signal to the input u = u + asin^cot), of 
small amplitude a and high frequency uj, and capturing an estimation of the gradient of the objective function 
by filtering. Three time scales (excitation signal, process, filtering) then operate in the overall closed loop 
system. Such a technique has been developed first in |49j for unknown functions ^(•) of Monod or Haldane 
type. Later on, an adaptive extremum-seeking scheme has been proposed to improve the performances of 
the transient response, for the Monod's case in [55] and the Haldane's one in [35]. In these approaches, the 
production rate r is assumed to be measured on-line, and explicit expressions of the unknown function /^(•) 
are required. An alternative extremum-seeking scheme combined with a neural network has been proposed 
in [20] to get out of these requirements. Nevertheless, all the mentioned approaches require the on-line 
measurement of the whole state (s, 6), excepted in [S^ where b only is measured but the result is dedicated 
to the Monod's kinetics. In the present paper, we consider that the single on-line measurement 

y = s (2) 

is available. One may consider coupling the former extremum-seeking controllers with state observers [171 
ini dni HH ISHl ISSl [Ml [2S] , but there is here an intrinsic difficulty due to the fact that "good" observers, i.e. 
whose speed of convergence can be arbitrary tuned, require the exact expression of the dynamics, and that 
"fast" observers are often sensitive to measurement noise and input disturbances. 

In the present work, we present a rather different approach that does not require any knowledge on the 
growth function neither the measurement of the objective function, and that is based on a continuation 
method [T]. We do not consider any excitation signal. Instead, we design a dynamical extension of the 
chemostat model ([T]) in a slow- fast or "singularly perturbed" form [53], such that the attracting critical 
manifold is exactly the graph of the unknown function s i— /i(s). Our extremum seeking scheme is an 
extension of a former work ^11 SI] that gives such a method for the reconstruction of the graph of /i(-) with 
the single measurement of the substrate, without any a priori knowledge on /i (excepted to be a smooth 
function). In particular, this technique can cope with non-monotonic growth functions and allows to identify 
unstable states of the open-loop system. The paper is organized as follows. In Section [5] we first present 
our general methodology, and then show how to apply it on the chemostat model in Section |5j Section [6] is 
devoted to numerical demonstrations. 



2 Assumptions in the general framework 

We consider a single- input / single-output dynamical system 

X = f{x,u) , . 

y = Kx) 

{x(t) £ M", u{t) £ M, y{t) £ M) and an objective function 

z = 0(y) (4) 

to be optimized at steady state, i.e., we are looking for an output feedback controller that steers the state 
to an operating point {x*,u*) that fulfills 

f{x\u*)^0 

(t>{h{x)) is locally max. w.r.t. m € R at a;* 
where the functions / and are unknown or partially known. 

We follow the statement of the general nonlinear problem of extremum seeking given in 2^ , that we adapt 
here to an output feedback framework. 
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Assumption Al (Existence of stabilizing output feedback). There exists a smooth output feed- 
back u{t) = a{y{t),u,y) parametrized by the parameter pair {u,y) €U xy (a pair of reference input and 
output), such that 

a{y,u,y) = u, (6) 

and the closed loop system 

X = f{x,a{h(x),u,y)) (7) 

admits a unique equilibrium X(.q{u,y), that is locally asymptotically stable for any (m, y) £U xy. 

For example, for the chemostat (H]), ^ the output feedback could be of the form u{t) — u + Gi{s ~ s{t)) 
with a sufficiently large Gi \42\ . 

Under Assumption Al, we consider the extremum-seeking problem for the closed- loop dynamics ([7]) as 
if the pair (m, y) was a new control, and look for pairs such that h{xeqiu, y)) = y- For this purpose, we shall 
consider the output-input characteristics defined as follows. 

Assumption A2 (Output-input characteristics). There exists a smooth function ip : y i-^ U such 
that 

y = h{xeq{u,y)) ^ u = ^:{y) (8) 

for any y € y. 

This assumption means that for each parameter y € y, there exists an unique input u £ U that is a zero 
of the function u h{xeq{u,y)) — y. For instance, in the chemostat model, positive equilibriums have to 
fulfill a{seq,u, s) = fi{seq), and then having Seq = s, along with ([B]), amounts to write u = ii{y) . 

Assumption A3 (Existence of local msLximum). There exists y* y such that the function (/){■) 
has a local strict maximum at y* . 

Finally, we shall assume that the unknown objective function (/>(•) possesses a known structure. 

Assumption A4 (Structure of objective function). There exists a smooth function (/? : i-> K 
such that 

^{y)=cp{y,^P{y)) yyey (9) 

Assumptions A2 and A4 imply that the extremum seeking problem amouts to optimize a known function 
of y and u at steady state. Thus, the problem consists then in finding an output feedback strategy that 
drives the state of the system to x* such that h{x*) — y* using the knowledge of the feedback law a and the 
objective only. The following two sections present two approaches. 



3 Continuous- time adaptation 

The continuous-time approach first constructs a differential equation for u that achieves u — ip{y) (and, thus, 
y — h{x)) asymptotically, and then applies a gradient search (again through a differential equation) along the 
curve of (u, y) given implicitly by y = h{x{u, y)). Note that for feedback laws of the form u{t) = u+Gly—y^t)] 
this implicit curve corresponds to the curve of open-loop equilibrium outputs. 



3.1 Step 1: an adaptive scheme for u 

We look for a dynamics 



with u, y) = 0, such that 



u = P{y,u,y) 
xeq{ip{y),y) 



(10) 

(11) 
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is a locally asymptotically stable equilibrium of the coupled dynamics 

i = f{x,a{h{x),{u,y))) 
u = P{h{x),u,y) 

Note that Assumptions Al and A2 ensures that E{y) is an equilibrium of (IT^ for any y G y, and that one 
has 

lim y{t) = y . (13) 

r— >+oo 

This is a classical adaptive output control problem, for which several techniques are available in the literature 
[5J [21] ■ We simply require the convergence 

lim u{t) = ip{y) (14) 

to be uniform w.r.t. y € y. 

3.2 Step 2: a continuation method for y 

Step 1 with Assumption A4 provides an approximation of the unknown objective function 0(-) at a fixed y: 

Hy) - <i>{y) = "piviu) . (15) 

The continuation consists in having y evolving slowly (step-wise or continuously) with a steepest descent 
based on an estimation of the gradient of the function (/). A continuous adaptation can be written as 



y 



di(p{y,u) + d2(p{y,u)ip'{y,u) (16) 



where e > is small compared to the time scale of dynamics (1121) . and ip'^y^u) is an estimation of the 
derivative of the function ip at y, for which u is an estimation oiip{y). Several gradient estimations techniques 
are available in the literature [301 ISIl EH UHl IS] • In the scalar case, an estimation of the sign of the derivative 
of the function is enough to choose the steepest descent. We can use for instance a dynamics with delay: 

= S{y,u,y{t-T),u{t-T)) 

= sgii[{(p{y,u) - (p{y{t-T),u{t-T)){y -y{t-T))] . 

Then, the overall dynamics 

i = f{x,a{h{x),u,y)) 

U = P{h{x),u,y) (18) 
y = e6{y,u,y{t-T),u{t-T)) 

is slow-fast where u — ip{y) is the (attracting) critical manifold. 

Finally, the output feedback strategy that we consider takes the following form 

u = a{y, u, y) with | " " ^j.^^ ""l ^1,^ , (19) 
1 y = t5(y,u,y(t-T),u{t-T)) ^ > 

Let us underline that the output y is naturally filtered by the dynamics to obtain the estimation y* as 

r= lim y{t) (20) 

that provides a robustness w.r.t. to measurement noise. Of course, the price to pay is to have a slow 
convergence due to small e. 
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4 Extremum-seeking using numerical optimization algorithms 



We compare the continuous-time adaptation with classical numerical optimization algorithms that are diffi- 
cult to express as dynamical systems such as golden-section iteration (which is applicable for single-parameter 
problems) and Newton iterations. These algorithms do not require any knowledge of the underlying system 
apart from the ability to evaluate the asymptotic output limt_j.oo y{t) for any given admissible inputs {u, y). 
Once the stabilizing feedback law uit) = a{y{t),u,y) is implemented, the problem is reduced to a classical 
numerical optimization problem: <j)[h[xeq{u, y))) — >■ max. In fact, if the feedback law is linear of the type 

u{t) = a{y{t), u, y)=u + G[y - y{t)] (21) 

this is an optimization problem in the single parameter v = u + Gy: 

(j){h{xeq{u + Gy))) F{v) ~> ma.x w.r.t v. (22) 

We test a combination of two classical optimization algorithms. We start with an initial admissible interval 
[viow , Vup] ■ Then we use golden section search to narrow the initial interval down to a given fraction of its 
original size and then apply a Newton iteration to solve F'{v) = in the remaining interval. Neither of these 
two methods requires knowledge about the structure of F, only the ability to evaluate F and its derivatives 
at any desired point (where derivatives can be approximated by finite differences). For the evaluation of F 
and its derivatives in a given point vq we use a naive act-and-wait approach [23] : one sets the parameter 
V — vq, such that the feedback control law is u{t) = vq ~ Gy{t), waits until the transients have settled 
(such that X = xq := Xeq{vo)), and then reads oS y = h{xo) and evaluates (j>{y) — ip{y,VQ — Gy). For the 
golden section search we iterate from an initial interval [wi,zj3] by setting V2 = vi + {vs — vi){3 — V5)/2 
and evaluating Fi — F{vi), F2 — F(v2) and F3 — F(v^) to obtain an initial triplet {vi,V2,vz) and the 
corresponding objective function values (i^i,F2,F3). Then we proceed with the standard golden search to 
iteratively obtain new triplets (vi,?;2,W3) until V3 — vi < tol. Correspondingly, for the Newton iteration 
with an initial guess Void we evaluate F at the three points {void — h, VoM, Void + h) (where his a. small finite- 
difference deviation) and choose as the new point Wnow the maximum point of the interpolating parabola. The 
golden section search is known to converge globally for unimodal functions with rate (-^5— l)/2, whereas the 
Newton iteration converges quadratically close to a local maximum. In the presence of random disturbances 
in the evaluation of F we gradually increase the accuracy of the evaluation of F during the iteration by using 
the mean of y over a longer time interval after the transients have settled. Note that, when using act-and- 
wait and a discrete-time optimization algorithm one has removed two of the two slower time scales present 
in [1], replacing them with the algorithmic iteration (the Newton iteration converges super-exponcntially) . 



5 Extremum seeking for the chemostat 

Reference [42j has shown that there exists an output feedback that satisfies the assumptions Al, A2 and A4 
for the model ([T]), with a simple saturated proportionate controller. Define first the saturation function as 
follows. 

e if Ce [Dmi„,^max], (23) 

^min if ^ <^ ^min? 



sat[D„,i„^n„„,](0 



For this saturation function [42 proved the following proposition. 

Proposition 1 Assume that the reference values {s,D) are in [smimSin) x [i'min, i'max]; where the bounds 
satisfy 

Dmin < m(s) for all s G [s 

-Dmax > /or a?/ S e [0, Sin]. 

Then the output feedback 

u{y, D, s) = sa<[z,.„,„,o.„.] {D ^ Gi{y ^ s)) (25) 
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with a gain Gi satisfying 

Gi > max ( max -^i'{s) , R^A^] (26) 

V6e[0,Si„] Sin - s / 

guarantees that the closed-loop dynamics ([^51) has a stable equilibrium (seq,6eq) G [0,Sin) x (0,oo), which 
attracts all initial conditions (s(0),6(0)) € [0, Sjn) x (0,cx)). 

See [42] for a proof. 

From system ([T]), one has at steady state 

/x(Seq) =5-Gi(Seq-s), (27) 

and consequently one has 

Soq = S ^ D ^ l^{Scq), (28) 

that is, Assumption A2 is fulfilled. 

Note that in system ([Ij at steady state one has 6 = Sin — s. Thus, the objective function 

(/)(s) =/i(s)(sin-s), (29) 

which is non-negative such that 0(0) = ipisin) — 0, is equivalent to the production rate ^{s)b. For </> the 
assumptions A3 and A4 are fulfilled. 

5.1 Implementation of continuous-time adaptation 

For the continuous-time adaptation we can choose the same adaptive scheme for u for the step 1 described 
in Section [01 as was chosen in [55] for the identification of /i. 

Proposition 2 For any s G (0,Si„) and numbers I?min, ^max such that < I?min < A*(s) < ^max; the 
dynamical feedback law 

u{y,D,s) = sat[D„^^^^^D,,-,^^]{D - Gi{y - s)) ^ 

D = -G2(2;-S)(^- Ani„)(Anax-5) ^ ' 

exponentially stabilizes the system ([T]) locally about (s, Sm — s), for any positive constants (Gi, G2) such that 
Gi > — /i'(s). Furthermore one has 



hm D{t)=ti{s) 

t—^+oo 



See [H] for the proof. 



Our continuous-time extremum seeking output feedback for the chcmostat model ([T]) is thus given by the 
following equations 

u{y, D, s) = sat[£,^.„ {D - Gi{y - s)) (31) 

b= -G2(y-S)(5- Anin)(i?max-5) (32) 

i| = sgn [(^(s,„ - s) - ^(i - r)(s„ - s(t - t)))(s - s(t - r))] (33) 
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5.2 Implementation of discrete-time optimization 



The discrete-time optimization uses only the saturated feedback, combining the parameters D and S into a 
single parameter v — D + Gis, such that 



u(y, v) = sat[ 



Giy) 



(34) 



If the initial bracketing interval [vi, V2] for v is chosen such the line £{v) = {{D, s) : D + Gis — v} intersects 
the domain of admissible reference values [smin, Sin) x [-Dmin, ^max] for all v S [wi, ■112], Proposition [T] ensures 
the existence of a unique equilibrium Soq(w) for all w G [vi, V2\- Thus, we can proceed with the golden section 
search for the objective functional 



F(v) 



u{Scqiv),v)[si, 



(35) 



where u, as defined in (|34p . is the asymptotic value of the input u. Once, the golden section search has shrunk 
the bracketing interval to a certain size, we apply a single step of the Newton iteration for the objective 
functional F. 



6 Numerical simulations 

The output feedback law (PT|) - ([55)) has been simulated for a Haldanc function 



K + s + s'^/Ki 



(36) 



and the following values of the parameters: 







Sin 


^min ^max 


1.0 1.0 


0.1 


1.0 


0.0 1.0 



Note however, that in our simulations we make no assumptions about the type of fi apart from those required 
in the assumptions of propositions [T] and [2] (leading to the validity of assumptions A1-A4). For the control 
laws in (|5T |l - (P5)) . the following parameters have been chosen: 



Gi G2 



1.0 1.0 le- 



To simulate measurement noise, we have considered that the output is corrupted in the following way 

yit)^sm + m) (37) 
where b{-) is a random square signal of frequency uj and amplitude a, whose values are given below. 



0.2 0.05 



6.1 Continuous-time adaptation 

Figure [1] shows how the trajectory of the variables (s, D) converges to the graph of the function i nthe 
(s, D)-plane initially up to disturbances due to the measurement noise. This convergence toward /i(-) is 
effected by the control laws ((3T|) and p2|) on the fast timescale (for t ^ 100, not visible in the time series 
depicted in Figure [2]). Once the feedback laws pi|) and ([32l) have achieved quasi-equilibrium for fixed s the 
dynamics of (1551) maximizes the functional (p on the slow timescale. Timescale separation e = 10~^ is chosen 
for our simulations. In Figures [T] and |3] one can also see that the adaptation of s via (j33p forces the state 
(s, b) to converge near the equilibrium (s*, b*) (indicated as intersection between dotted line and graph /z(-), 
where s* is maximizing the function </>(•). We note that the adaptation has converged already after t — 3500. 

Remark If we choose the timescales not sufficiently different (say, e = 10~^ in our example) we observe 
that the chemostat ([T} with feedback laws ([5T |) - ([551) becomes dynamically unstable in in its fixed point 

{s,b,^s,D) = is\b\s\^i{s*)). 
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Figure 2: Output of JT]) with feedback laws ([3T1) ~ ([33| . disturbed by ([BT)) . over time (same run as Figured]). 



6.2 Discrete-time algorithms 

For the discrete-time optimization algorithm we chose an initial bracketing interval [wi,t'2] = [0.04, 1] (gain 
Gi = 1 in feedback law (|34l) identical to the continuous-time adaptation). During the golden section search 
(until t « 700 in Figure ID^a)) we checked every ijnc = 25 time units if the transients have settled (criterion is 
a decrease of the standard deviation compared to the previous time interval) . Once the output is accepted as 
stationary we used the mean of output and input over the last time interval to calculate the objective F. A 
larger tine inceases the accuracy of the evaluation of F in the presence of disturbances due to averaging. We 
observe that after 7 evaluations of F the golden section search reduces the bracketing interval to length 0.2 
and the estimated optimal value of s is already close to s* {t « 700 in Figure HJ^a)). The Newton step then 
only confirms the optimality of the final estimate of s* (up to disturbance level). For the Newton iteration 
we chose ti^c larger {ti^c — 100) to increase the accuracy of F evaluations and use the points (v — h,v,v + h) 
with h = 0.05 to evaluate the approximating parabola for F in v. Figure ID^b) shows the trace (grey dots) 
of outputs s and inputs u (equaling dilution rate D) in the (s, Z3)-plane. The gaps in the trace are due to 
finite sampling time. The full circles show the points where the transients were accepted as settled. 

The speed and accuracy of the discrete-time optimization is limited by the choice of iinc, which in turn 
has to be determined sufficiently large to enable an averaging effect for the disturbances. 
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Figure 3: Trajectory of same run as Figures [1] and [5] in the (s, b) plane {{{s,b) : s + b — Sin} is invariant). 

7 Conclusion 

We have proposed two extremum-seeking schemes: one approach uses dynamical output feedback that is 
based on a continuation method. Its closed-loop dynamics possesses two time scales: the timescale of the 
original system with stabilizing feedback and a slow timescale for the adaptation scheme. As an alternative we 
test classical optimization algorithms, applied to the steady-state outputs of the feedback-stabilized system, 
which have a continuous timescale (which is exponentially convergent) and a discrete timescale (which is 
exponentially or superexponentially convergent). This is in contrast to extremum seeking techniques that use 
dither signals, which need three timescales for nonlinear systems (the time scale of the stabilized system is 
the fastest). We illustrate our approaches on the chemostat model for unknown growth function and single 
measurement. The numerical simulations suggest that the method is practically viable and robust with 
respect to measurement disturbances. Extensions and a detailed convergence analysis for general systems 
and the continuous-time adaptation will be the subject of future work. Conditions for the convergence of 
the discrete-time optimization can be deduced directly from the corresponding statements for the original 
numerical algorithms. 
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